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The variable range hopping theory, as formulated for exponentially localized impurity states, 
does not necessarily apply in the case of graphene with covalently attached impurities. We analyze 
the localization of impurity states in graphene using the nearest neighbor tight-binding model of 
an adatom-graphene system with Green's function perturbation methods. The amplitude of the 
impurity state wave function is determined to decay as a power law with exponents depending 
on sublattice, direction, and the impurity species. We revisit the variable range hopping theory 
in view of this result and find that the conductivity depends as a power law of the temperature 
with an exponent related to the localization of the wave function. We show that this temperature 
dependence is in agreement with available experimental results. 



Chemical functionalization of graphene has been pro- 
posed as one of the most promising methods to modify 
its transport properties. The attachment of covalently 
bonded atoms or functional groups opens the possibility 
to design devices [HE], open band gaps 014], and control 
the excellent transport properties of its massless Dirac 
fermions [5]. Because graphene is essentially a two di- 
mensional electronic system, even small amounts of this 
functional groups can produce radical changes into its 
transport properties. The temperature dependence of 
the conductivity of chemically functionalized graphene 
shows a peculiar behaviour that strongly suggests the 
importance of disorder. This phenomenon is very gen- 
eral and has been observed with hydrogen [S] [7] , fluorine 
ES M i oxygen [TMT2] and metals [T3] ■ In order to relate 
the experimental measurements with microscopic prop- 
erties, this anomalous temperature dependence has been 
analyzed with a variable range hopping (VRH) theory as 
formulated for semiconductors with exponentially local- 
ized impurity states [TH [T5] . As a consequence, the esti- 
mated localization length does not seem to correlate with 
any reasonable characteristic length in these systems. In 
the case of dilute fluorinated graphene, the localization 
length is obtained to be 56 nm [9] , while it is estimated to 
be 90 nm in lightly silver-coated graphene [13] ■ The VRH 
theory, as originally formulated by Mott, is not applicable 
when the localization of the impurity states is not expo- 
nential. For the case in hand, it has been extensively 
discussed that the impurities form resonant states and 
the low temperature conductivity as a function of carrier 
concentration has been theoretically determined in good 
agreement with experimental results [T6HI9] . Here, we 
determine the power law decay of these impurity states 
and show that the exponent depends on the resonant 
energy and approaches asymptotically the case of vacan- 
cies [2D1 HJJ . The exponent is also anisotropic displaying 
strong dependence on the sublattice and on the direction. 
We use our findings to reformulate the VRH theory under 
the assumption that there is a regime of temperature and 
density where the jump between these impurity states 
is incoherent and determine a general behavior that ex- 



plains the measured temperature dependence of the con- 
ductivity in these systems. The use of this reformulated 
theory provides a method to determine the localization 
characteristics of the impurity states in graphene. 

We start by studying the system of a single impurity 
atom on graphene with a tight-binding Hamiltonian of 
a localized p z -orbital basis set \i) for the it band of the 
pristine graphene and |ad) for the adatom, 
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where t is the hopping energy between nearest neigh- 
bor carbon atoms (ss 2.8 eV). e a d is the site energy of 
the adatom and V&d is the hopping energy between the 
adatom and the carbon atom to which it is attached (|0)). 
This model is often used to study the adatom-graphene 
system [TBI E21 E3] . Treating H' as a perturbation, the 
T matrix is given by 
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where G is the Green's function of H , G%° = (0| G |0) 



and Gg' 
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The per- 



turbed eigenstate in the band continuum is given by the 
Lippman-Schwinger equation 
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For a certain energy E — e r that satisfies the resonance 
condition 



Re[l-V^(e r )G^(e r )} = 0, 
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FIG. 1. Amplitude of the resonance state at different resonance energies, (a) e r — t/300. The amplitudes on the A sublattices 
vanish, (b) e r — t/6. The insets show the amplitudes represented by the radii of circles on the graphene honeycomb lattice. 
In both plots a circle symbol represents an amplitude on a B-site, which is in the direction indicated by a line with the same 
color in the insets. A square symbol represents an amplitude on an A-site, which has a circle drawn with the same color in the 
inset. The center dot (blue) is where the adatom is attached, a is the nearest neighbor carbon-carbon distance. 



the second term in Eq. ^ would be significantly en- 
hanced and \ipo(E)) can be ignored near the impurity 
site, which gives 



(i | i){er)) cx (t | G {e r ) | 0) . 
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This resonance state will have larger amplitude near the 
adatom but it never decays to zero for large distance be- 
cause of the contribution of the Bloch function \ipo(E)). 
Also, note that this expression only depends on the reso- 
nance energy e r , which means the values of e a d and in 
the Hamiltonian only contribute to determine e r through 
Eq. Q and we can study the effects of different adatom 
species by varying e r . 

The decay of the wave function of the impurity state 
can be studied by investigating the lattice Green's func- 
tions (GFs) (i | Go(e r ) | 0), which are determined mostly 
from contributions from the two Dirac points (K, K') 
in the Brillouin zone when they are evaluated at ener- 
gies near zero. Integrating around the two Dirac points 
rather than the whole BZ and assuming a completely lin- 
ear band, the GFs are calculated and given in terms of 
Hankel functions j2"H-j2~o] (labeled A for sites in the same 
sublattice as the impurity and B for sites in the opposite 
sublattice) 
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where A c is the area of a unit cell in graphene and vf 
is the Fermi velocity. The amplitude of the Hankel func- 
tions Hq 1 ^ and decay isotropically, but the GFs also 



depend on the prefactors 
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P = e^ r + e lK ' r , (11) 

where r = tan -1 (r y /r x ) when the x-axis is taken to be 
along K' — K. The form of the argument of the Hankel 
function makes two types of approximations feasible. For 
an impurity with a small e a d or a large T4d (e.g. a va- 
cancy), the resonance energy e r solved from Eq. [6] will be 
small, which means we can do small argument expansion 
to the Hankel function. This gives a resonance state that 
has zero amplitude on the A sublattice sites and decays 
as r _1 for the B-sites [5DJI2S]. On the other hand, when 
e r is not vanishingly small, we are more interested in the 
long range decaying behavior and it is necessary to do 
large argument expansion, which gives 



v F J 



(2v F \ 
\wEr J 



1 



[Av 2 - l) v F 
8Er 



(12) 

Given enough distance, both the A-site and the B-site 
amplitude will fall off primarily as r~ 0,5 . 

The decay behavior is further elucidated by evaluating 
the GFs directly. The method used here to obtain the 
GFs for the honeycomb lattice follows a calculation for 
square lattice [57] , in which the lattice GFs are calculated 
from larger to smaller distances from the impurity. This 
is done to avoid a diverging term which originates from 
numerical instabilities and also satisfies the GF equation 
of motion [281129] , 

The calculated GFs (amplitudes of the resonance 
state) are plotted in Fig. [l] where in the insets the am- 
plitude of the GF on a certain site is represented by the 
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FIG. 2. The characteristic decay exponents of the two sub- 
lattice sites vs. the resonance energy. For the A-sites, the 
exponent is taken from the sites that forms a triangular lat- 
tice with the impurity site, similar to the top line in Fig.[TJb). 
For the B-sites, the exponent is taken from sites in the arm- 
chair direction, similar to the green line in Fig. [T] 



radius of the circle that is drawn on that site, and the 
results mostly confirm the approximations. Firstly, the 
amplitude of the resonance state wave function depend 
on the resonance energy e r and the energy dependent be- 
haviors of the two sublattice sites are drastically different. 
This is clear from Fig. Qa) (e r = t/300) and Fig. Qb) 
(e r = t/6). At low resonance energy, the resonance state 
is almost exclusively on the B sublattice sites. The am- 
plitudes on the A-sites increase quickly while the B-sites 
stay relatively the same with increasing e r and the two 
sublattice sites have comparable amplitudes at e r = t/6. 
Secondly, the wave function amplitude decays with power 
law 
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although the exponent s depends on resonance energy, 
sublattice and direction. For the A-sites (Fig.JlJb)), the 
sites that form a triangular lattice with the impurity site 
(marked with red) have larger amplitudes than the other 
sites (marked with black) when they are approximately 
the same distance to the impurity. This can be explained 
with the prefactor /3 in Eq. which evaluates to 2 for 
the former group of sites and -1 for the latter. The two 
groups give almost perfect linear fits in a log-log plot 
close to the impurity with essentially the same decay ex- 
ponent s, while the second group of sites deviate from 
the linear fit at larger distance. For the B-sites, the de- 
cay is anisotropical and power laws can be seen in many 
directions when e r is small (Fig. [lja)). This behavior 
has been obtained with approximations to the GFs pre- 
viously by Nanda et al. [2B], and our calculation confirms 
their result. At higher e r , the decay in the B sites in the 
armchair direction is the slowest (thus contributes the 



most when calculating overlap) and still obeys very good 
power laws, while the other data sets start to deviate. 
Deviations from power laws in both the A sites and the 
B sites are not explained in the approximated GFs and 
they happen at a smaller distance for larger e r . There- 
fore, it is possibly a result of non-zero energy and contri- 
butions from k-points other than the two Dirac points. 
The decay exponents for the A-sites and the armchair 
direction of the B-sites are plotted in Fig. [2] with several 
resonance energies. As predicted by the approximations, 
the B-sites decay exponent is 1 at zero resonance en- 
ergy [501111] and both exponents approach 0.5 with large 
resonance energy. Finally, we would like to note that 
the amplitude of the wave function calculated with GFs 
agrees with evaluation of the GFs with elliptic integrals 
[28] . the results in Ref . |2"31 for vacancies (e r = 0), as well 
as the amplitude of the eigenfunction near the impurity 
(r < 10a) obtained from direct diagonalization of the 
Hamiltonian with periodic boundary conditions. 

The power-law decay and the Bloch-wave behavior at 
large distance both indicate that the impurity state in 
graphene is not nearly as localized as a typical midgap 
state in a semiconductor, which decays exponentially. 
The resonance state here are not normalizable and it 
would be difficult to define a localization length. In view 
of the simplicity and the extensive usage of the VRH the- 
ory, it is necessary to investigate the impact of a power- 
law-decaying impurity state on the hopping conductivity 
result. 

Similar to VRH, we will assume that the density of 
impurities is low enough so that the average distance be- 
tween impurities is longer than the phase coherent length 
and coherent scattering from multiple centers can be ig- 
nored. The derivation of VRH is outlined in Ref. [F4l and 
loTJl By simply replacing the overlap with Eq. (13 1, the 



hopping probability between two impurity sites i and j 
can be written as 
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where and are the distance and energy difference 
between the two states, respectively. 7^ is a prefac- 
tor that comes from the coupling between electrons and 
phonons. It depends on ey and with power laws and 



is thus ignored in the original VRH 



Assuming a 



smooth density of states near the Fermi level (<7(ep)), £j 
and Tij are related by 



-l/d 



(15) 



where d is the dimension of the system. Taking the oc- 
cupation of the two states into account, we obtain the 
conductance between two impurities as (301 13 1] 
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Conventional VRH 


Power Law VRH 


Data Set 


r 2 




r 2 




Hydrogenated, V g = OV, Ref. M 


0.9857 


284 


QQ30 


0.714 


Hydrogenated, V g = OV, Ref. [7J 


0.9884 


280 


0.9975 


0.524 


Hydrogenated, V„ = 3V, Ref. El 


0.9819 


187 


0.9981 


0.459 


Hydrogenated, V 9 = 9V, Ref. [7J 
Fluorinated, V g = OV, Ref. [S] 


0.9621 
0.9843 


107 

450 


0.9928 
0.9730 


0.382 
0.826 


Fluorinated, V g = OV, Ref. [9] 
Fluorinated, n = 0.7 x 10 12 cm" 2 , Ref.H 


0.9990 


270 


0.9761 


0.664 


0.9985 


130 


0.9843 


0.550 


Fluorinated, n = 1.4 X 10 12 cm" 2 , Ref.H 


0.9889 


33 


0.9955 


0.335 


Fluorinated, n = 2.5 x 10 12 cm -2 , Ref.H 


0.9599 


5 


0.9929 


0.190 



TABLE I. The fitting details of available data sets with both the conventional VRH and Eq. | |17| . V g is the experimental gate 
voltage relative to the charge neutrality point (the voltage where the sample exhibits the highest resistance), n is the charge 
density calculated with the gate voltage and sample specifics. To is the characteristic temperature in the conventional VRH 
In a oc — (Tb/T) 1//3 and r 2 is the squared correlation coefficient for the linear fits. 



where e is the charge of an electron. The conductivity of 
the bulk is assumed to be proportional to aij between the 
most conductive pair of impurities, which corresponds 
to the maximum of cry when varying or tij. This 
leads to a power law temperature dependence for the 
conductivity, 



a oc T v , with r\ 



2s 



(17) 



where s' originates from the prefactor jij/T. For hydro- 
genic states, it is estimated to be {v — 2)/(d+ 1), where 
v is the critical exponent for the size of the percolating 
cluster [30]. In 2D, v = 1.34 [32] and 



-0.22 



(18) 



The analysis that gives this exponent cannot be easily 
generalized, but we can still assume s' to be a constant 
that does not depend on the details of the impurities. 



Eq. (17) can be used to fit existing experimental data 



of the systems of hydrogen adatoms [6j [7] and fluorine 
adatoms on graphene [5] and the extracted parameters 
are shown in Table [I] along with fits of the original VRH. 
The fits are done to all the data points presented in the 
references. Assuming Eq. 18 the exponents extracted 



from the fittings are within the reasonable range that is 
expected from this theory. In both experiments where the 
effect of gate voltage is studied, the exponent decreases 
when the gate voltage moves away from the charge neu- 
trality point, effectively shifting the Fermi energy so that 
on average impurity states with higher resonance energies 
participate in the conduction. This is consistent with the 
behavior of the B-sites decay exponent (Fig. [2]) while the 
amplitudes on the A-sites are too small to be relevant in 
the range of the experimental voltage. 



The conventional VRH and Eq. (17) have very similar 



curvatures in this temperature range and all the data can 
fit both equations fairly well. Comparing the correlation 
coefficients, the data for fluorinated graphene at low gate 



voltage fit the original VRH better, while the power law 
dependence is more pronounced in all the other data sets. 
The currently available data cannot convincingly exclude 
either equation as the conduction mechanism, especially 
considering the VRH requires low temperature but all the 
data go up to room temperature. We expect more contin- 
uous and accurate experimental data at low temperature 
to prove our proposal of a power law temperature depen- 
dence. 

In summary, we have shown that the impurity state 
in graphene is a resonance state in the band continuum 
and it is localized only as power law functions with ex- 
ponents generally below 1. This means that the VRH 
theory which assumes exponential localization is not di- 
rectly applicable to disordered graphene. Replacing the 
overlap term in VRH, a theory for the temperature de- 
pendence of conductivity is derived which fits the existing 
experimental data. However, since the states are largely 
delocalized, the hopping picture of conduction may not 
be the most appropriate approach to model the transport 
properties of these systems. Further investigation into 
this problem is needed to develop a theory that includes 
both the impurity states and the extended unperturbed 
states. 
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